#include <iostream>
#include <fstream>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <functional>
#include <cmath>
#include <vector>
#include <utility>
#include <map>
#include "Types.hpp"
#include "Vertex.hpp"
#include "Colsamm.h"
#include "sMatrix.hpp"
#include <unordered_map>
#include <valarray>
#include "VTKWriter.hpp"

#define PI 3.14159265
using namespace ::_COLSAMM_;
real delta;

sMatrix readVertexFile(const std::string& fileName,std::map<size_t,Vertex>& vertex,std::vector<Face>& faces, int level_global, int level_local);
std::vector<Face> refineMesh(int level, std::vector<Face> &faces,std::map<size_t,Vertex> &vertex );
std::vector<Face> LrefineMesh(int level, std::vector<Face> &faces,std::map<size_t,Vertex> &vertex );

int main(int argc, char *argv[])
{
    if (argc !=3 && argc!=4)
    {
        std::cout<<"Please enter the right number of parameters\n ./waveguide delta epsilon \nProgram exiting...\n"; // Valid even if the Debug mode is Off.
        return 1;
    }
    int level_local=0;
    //Parse the parameters
    real epsilon = atof(argv[1]);
    int level_global = atof(argv[2]);
    std::cout << "GLOBAL REFINEMENT REQUIRED UPTO LEVEL:\t" << level_global << "\n";
    if (argc ==4)
    {
        level_local = atof(argv[3]);
        std::cout << "LOCAL REFINEMENT REQUIRED UPTO LEVEL:\t" << level_local << "\n";
    }
    else std::cout << "No Local Refinement Requested" << std::endl;
    //Map of index to Vertex
    std::map<size_t,Vertex> vertex;
    //Vector of faces
    std::vector<Face> faces;
    //read the file
    sMatrix sparseStiffnes = readVertexFile("L_unit.txt",vertex,faces, level_global, level_local);

    ELEMENTS::Triangle my_element;
    std::vector< std::vector< double> > local_stiffness;
    std::vector<  double> corners(6,0.0);

    for(auto iter = faces.begin(); iter !=faces.end(); ++iter )
    {
        size_t vertexId[3];
        for (size_t i = 0; i<N_FACE_VERTEX; ++i)
        {
            vertexId[i] = (*iter).vertexIDs_[i];
            corners[2*i]= vertex.find(vertexId[i])->second.getCoordinate(0);
            corners[2*i+1]= vertex.find(vertexId[i])->second.getCoordinate(1);
        }
        my_element(corners);

        local_stiffness = my_element.integrate(grad(v_())*grad(w_()));

        // Local Stiffness Matrix is linked to the vertexId[i,j,k] . Need to assemble the Global Matrix through it
        for (size_t i=0; i<3; ++i)
        {
            for (size_t j=0; j<3; ++j)
            {
                sparseStiffnes(vertexId[i],vertexId[j])+=local_stiffness[i][j];
  //             std::cout << vertexId[i] << " " << vertexId[j] << "value "<<  local_stiffness[i][j] << "\n" ;
            }
   //         std::cout << std::endl;
        }
 //       std::cout << std::endl;

    }
    std::ofstream file_outLocal("stiffness.txt");
    printMat(file_outLocal,sparseStiffnes);
    file_outLocal.close();

    //  for(auto f : faces)
    //   std::cout << f.vertexIDs_[0] <<" "<< f.vertexIDs_[1] <<" "<<f.vertexIDs_[2]<<"\n";
   // for(auto f : vertex)
   // std::cout << f.second.getCoordinate(0) <<" "<< f.second.getCoordinate(1)<<"\n";


    std::vector<real> u= GSSolver(sparseStiffnes,  epsilon,1000, vertex);
    VTKWriter VTK("seminar_result");
    VTK.WriteData(vertex, u);

    // PLOTTING The Analytical Solution vs. the Computed Solution
    std::ofstream file_solution("comparison.txt");
    file_solution << "X\t" <<"Y\t" << "U_analytical" << "U_Computed\n";
    real sum=0;
    for (size_t i =0; i< u.size() ; ++i)
    {
        file_solution << vertex[i].getCoordinate(0) << "\t" <<  vertex[i].getCoordinate(1) << "\t" << boundary1(vertex[i].getCoordinate(0), vertex[i].getCoordinate(1)) << "\t" << u[i] << "\n" ;
        sum+=(u[i]-boundary1(vertex[i].getCoordinate(0), vertex[i].getCoordinate(1)))*(u[i]-boundary1(vertex[i].getCoordinate(0), vertex[i].getCoordinate(1)));
    }

    std::cout<<sqrt(sum/u.size()) << std::endl;

    //for(auto f : u)
    //	std::cout << f <<"\n";
    return 0;
}




std::vector<std::string> tokenize(const char* str, const char* delim)
{ 
    std::vector<std::string> tokens;
    char* token = strtok (const_cast<char*>(str), const_cast<char*>(delim));
    while (token != NULL) {
        tokens.push_back(token);
        token = strtok (NULL, delim);
    }
    return tokens;
}

sMatrix readVertexFile(const std::string& fileName,std::map<size_t,Vertex>& vertex,std::vector<Face>& faces, int level_global, int level_local)
{
    std::ifstream infile (fileName.c_str(), std::ios::in);
    const int max_chars_per_line = 512;
    char *buffer(new char[max_chars_per_line]);
    std::vector<std::string> tokens;
    size_t n_tokens,N_vertex,N_faces;
    real coordinates[DIM];
    size_t id_vertex;
    if(!(infile.is_open() && infile.good()))
    {
        std::cerr << std::endl;
        std::cerr << "ERROR: Failed to open file: " << fileName << std::endl;
        std::cerr << std::endl;
        infile.close();
        return false;
    }
    infile.getline(buffer, max_chars_per_line);
    tokens = tokenize(buffer, " \t\r\n");
    N_vertex = atoi(tokens[0].c_str());
    infile.getline(buffer, max_chars_per_line);

    // Define the sMatrix based on the Number of Vertex - Constructor of the sMatrix could be modified

    for(size_t i=0;i<N_vertex;++i)
    {
        infile.getline(buffer, max_chars_per_line);
        tokens = tokenize(buffer, " \t\r\n");
        n_tokens=tokens.size();
        id_vertex=atoi(tokens[0].c_str());
        //Get coordinates
        if((n_tokens-1)!=DIM)
        {
            std::cerr << std::endl;
            std::cerr << "ERROR: Dimensions are not equivalent: " << DIM <<" != "<<n_tokens-1<< std::endl;
            std::cerr << std::endl;
            infile.close();
            return false;
        }
        for(size_t j=1;j<n_tokens;++j)
            coordinates[j-1]=atof(tokens[j].c_str());

        //creating new vertex
        Vertex vi(coordinates);

        std::map<size_t,Vertex>::const_iterator it=vertex.find(id_vertex),end_it=vertex.end();
        if(it==end_it)
            vertex.insert(std::pair<size_t,Vertex>(id_vertex,vi));
        else
        {
            std::cerr << std::endl;
            std::cerr << "WARNING: Two or more vertex have the same index: " << id_vertex << std::endl;
            std::cerr << std::endl;
        }
        //Setting the Sparse Matrix Outer Dimensions using the vertex Id's

    }
    //For the faces
    infile.getline(buffer, max_chars_per_line);
    tokens = tokenize(buffer, " \t\r\n");
    N_faces = atoi(tokens[0].c_str());
    infile.getline(buffer, max_chars_per_line);
    size_t face_vertexIds[N_FACE_VERTEX]; // Array of three Dimensions- One Face (with Id's)
    for(size_t i=0;i<N_faces;++i)
    {
        infile.getline(buffer, max_chars_per_line);
        tokens = tokenize(buffer, " \t\r\n");
        n_tokens=tokens.size();
        if(n_tokens!=N_FACE_VERTEX)
        {
            std::cerr << std::endl;
            std::cerr << "ERROR: The number of vertex in each face must be: " << N_FACE_VERTEX <<" not "<<n_tokens<< std::endl;
            std::cerr << std::endl;
            infile.close();
            return false;
        }
        for(size_t j=0;j<n_tokens;++j)
            face_vertexIds[j]=atoi(tokens[j].c_str());

        Face fi(face_vertexIds);
        faces.push_back(fi);
    }

    std::vector<Face> faces_refined= refineMesh(level_global,faces,vertex);
    faces.clear();
    faces= faces_refined;

    std::vector<Face> faces_Lrefined = LrefineMesh(level_local,faces,vertex);

    faces.clear();
    faces= faces_Lrefined;

    //New Neighbours need to be set with the Refined Mesh--Later?
    sMatrix sparseStiffnes(vertex.size());

    for (size_t it = 0; it<vertex.size(); ++it)
        sparseStiffnes.addRow(it);

    for (size_t it =0; it<faces.size(); ++it)
    {
        for(size_t c=0;c<n_tokens;++c)
        {
            for(size_t j=0;j<n_tokens;++j)
            {
                if(c!=j)
                {
                    vertex[faces[it].vertexIDs_[c]].addNeighborID(faces[it].vertexIDs_[j]);
                }
                sparseStiffnes.addElement(faces[it].vertexIDs_[c],faces[it].vertexIDs_[j],0);
            }
        }
    }

    // Set Boundary Flags for all the Vertices
    for (size_t it =1; it<vertex.size(); ++it)
    {
        vertex[it].setbdFlag();
    }

    vertex[0].setbdFlag(1); // Singular Point in this case is vertex0;
    delete [] buffer;

    infile.close();
    return sparseStiffnes;
}



std::vector<Face> refineMesh(int level, std::vector<Face> &faces,std::map<size_t,Vertex> &vertex )
{
    std::vector<Face> face_ref (faces);
    std::unordered_map<Edge,size_t, VertexSetHashFunction,SetEqual > edges;

    if (level >0)
    {
        std::cout << "Refining for level:\t" << level <<std::endl;
        face_ref.clear();
        int n_tokens =3;
        // Loop through The Faces in the vector<Face> faces
        for (size_t it =0 ; it< faces.size(); ++it)
        {
            // Splitting of the Face_Vertex Based on the Refinement Level
            size_t vertexIdCenter[3];
            for(int i=0;i<n_tokens-1;++i)
            {
                auto vertex0_iter = vertex.find(faces[it].vertexIDs_[i]);
                real coords[2];
                for(int j=i+1;j<n_tokens;++j)
                {
                    auto vertex1_iter = vertex.find(faces[it].vertexIDs_[j]);
                    Edge e1(std::pair<size_t,size_t>(vertex0_iter->first,vertex1_iter->first));
                    auto edge_iter = edges.find(e1);
                    if (edge_iter == edges.end())	{

                        Vertex v = vertex0_iter->second;
                        Vertex w = vertex1_iter->second;

                        coords[0]=(v.getCoordinate(0)+w.getCoordinate(0))*0.5;
                        coords[1]=(v.getCoordinate(1)+w.getCoordinate(1))*0.5;
                        Vertex center(coords);
                        vertexIdCenter[i+j-1] = vertex.size();
                        vertex.insert(std::pair<size_t,Vertex>(vertexIdCenter[i+j-1],center));
                        edges.insert(std::pair<Edge,size_t>(e1,vertexIdCenter[i+j-1]));
                    }
                    else
                        vertexIdCenter[i+j-1] = edge_iter->second;
                }
            }
            size_t refineFaceId[3];
            // Face 1
            refineFaceId[0] = faces[it].vertexIDs_[0];
            refineFaceId[1] = vertexIdCenter[0];
            refineFaceId[2] = vertexIdCenter[1];
            Face fi1(refineFaceId);
            face_ref.push_back(fi1);

            //Face2
            refineFaceId[0] = faces[it].vertexIDs_[1];
            refineFaceId[1] = vertexIdCenter[2];
            refineFaceId[2] = vertexIdCenter[0];
            Face fi2(refineFaceId);
            face_ref.push_back(fi2);

            //Face3
            refineFaceId[0] = faces[it].vertexIDs_[2];
            refineFaceId[1] = vertexIdCenter[1];
            refineFaceId[2] = vertexIdCenter[2];
            Face fi3(refineFaceId);
            face_ref.push_back(fi3);

            //Face4
            refineFaceId[0] = vertexIdCenter[0];
            refineFaceId[1] = vertexIdCenter[2];
            refineFaceId[2] = vertexIdCenter[1];
            Face fi4(refineFaceId);
            face_ref.push_back(fi4);
        }
        return refineMesh(level-1, face_ref,vertex);
    }
    else
    {
        std::cout << "Number of Faces are:\t" << face_ref.size()<< std::endl;
        std::cout << "Number of Vertices are:\t" << vertex.size()<< std::endl;

        return face_ref;
    }
}



std::vector<Face> LrefineMesh(int level, std::vector<Face> &faces,std::map<size_t,Vertex> &vertex )
{
    std::vector<Face> face_ref (faces);
    std::unordered_map<Edge,size_t, VertexSetHashFunction,SetEqual > edges;
    real sing_point[2]; sing_point[0]=0.0;sing_point[1]=0.0;
    Vertex Singular(sing_point);

    if (level >0)
    {
        std::cout << "Local Refining for level:\t" << level <<std::endl;
        face_ref.clear();
        int n_tokens =3;
        // Loop through The Faces in the vector<Face> faces
        for (size_t it =0 ; it< faces.size(); ++it)
        {

            // Check if this face has to be refined or not! Loop through the vertices of the Face
            bool ret=false;
            size_t p=0;
            while (!ret && p<3)
            {
                size_t vertexID = faces[it].vertexIDs_[p];
                //				auto vertex_iter = vertex.find(faces[it].vertexIDs_[p]);
                //				Vertex v = vertex_iter->second;
                //				if (v.isEqual(Singular)) { ret = true; break;} ++p;
                if (vertexID ==0) { ret =true; break;
                }
                ++p;
            }

            if (ret)
            {
                // Splitting of the Face_Vertex Based on the Refinement Level
                size_t vertexIdCenter[3];
                for(int i=0;i<n_tokens-1;++i)
                {
                    auto vertex0_iter = vertex.find(faces[it].vertexIDs_[i]);
                    real coords[2];
                    for(int j=i+1;j<n_tokens;++j)
                    {
                        auto vertex1_iter = vertex.find(faces[it].vertexIDs_[j]);
                        Edge e1(std::pair<size_t,size_t>(vertex0_iter->first,vertex1_iter->first));
                        auto edge_iter = edges.find(e1);
                        if (edge_iter == edges.end())	{

                            Vertex v = vertex0_iter->second;
                            Vertex w = vertex1_iter->second;

                            coords[0]=(v.getCoordinate(0)+w.getCoordinate(0))*0.5;
                            coords[1]=(v.getCoordinate(1)+w.getCoordinate(1))*0.5;
                            Vertex center(coords);
                            vertexIdCenter[i+j-1] = vertex.size();
                            vertex.insert(std::pair<size_t,Vertex>(vertexIdCenter[i+j-1],center));
                            edges.insert(std::pair<Edge,size_t>(e1,vertexIdCenter[i+j-1]));
                        }
                        else
                            vertexIdCenter[i+j-1] = edge_iter->second;
                    }
                }
                size_t refineFaceId[3];
                // Face 1
                refineFaceId[0] = faces[it].vertexIDs_[0];
                refineFaceId[1] = vertexIdCenter[0];
                refineFaceId[2] = vertexIdCenter[1];
                Face fi1(refineFaceId);
                face_ref.push_back(fi1);

                //Face2
                refineFaceId[0] = faces[it].vertexIDs_[1];
                refineFaceId[1] = vertexIdCenter[0];
                refineFaceId[2] = vertexIdCenter[2];
                Face fi2(refineFaceId);
                face_ref.push_back(fi2);

                //Face3
                refineFaceId[0] = faces[it].vertexIDs_[2];
                refineFaceId[1] = vertexIdCenter[1];
                refineFaceId[2] = vertexIdCenter[2];
                Face fi3(refineFaceId);
                face_ref.push_back(fi3);

                //Face4
                refineFaceId[0] = vertexIdCenter[0];
                refineFaceId[1] = vertexIdCenter[1];
                refineFaceId[2] = vertexIdCenter[2];
                Face fi4(refineFaceId);
                face_ref.push_back(fi4);

            }
            else  {
                Face fi_orig(faces[it]);
                face_ref.push_back(fi_orig) ;
            }


        }
        return LrefineMesh(level-1, face_ref,vertex);
    }
    else
    {
        std::cout << "Number of Faces are:\t" << face_ref.size()<< std::endl;
        std::cout << "Number of Vertices are:\t" << vertex.size()<< std::endl;

        return face_ref;
    }
}

